###########################################################
# COERCION PROVOCATION (2016) SURVEY EXPERIMENT CODE

rm(list = ls())

# print date and time of run
print("Starting date and time")
Sys.time()

# load packages and functions
R_folder <- "INSERT HERE"
source(paste0(R_folder, "survey_functions.R"))

# load data
data_folder <- "INSERT HERE"
load(paste0(data_folder, "coercion_provocation.RData"))

# images directory
images_directory <- "INSERT HERE"

# grayscale for the graphics or not
#print_colormodel <- "srgb" # in color
print_colormodel <- "grey"

# Function to make data frames for coefficient plots:
make.df <- function(var,comp,n,model,group){
  data.frame(
    var = as.character(var),
    comp = as.character(comp),
    n = n,
    coeff=model[2,1],
    se=model[2,2],
    ci90 = model[2,2]*(-qnorm((1-0.9)/2)), # 90% and 95% confidence intervals are calculated based on the normal distribution
    ci95 = model[2,2]*(-qnorm((1-0.95)/2)),
    pval = paste("p",ifelse(model[2,4] < 0.001,
                            "< 0.001",
                            paste("=", round(model[2,4], 3), sep=" "))),
    group = as.character(group)) 
  #one-sided p-value
}

# Function to estimate robust standard errors
coef_r <- function(lm_object, v_c_matrix = "HC1",
                   var = NULL, comp = NULL, group = NULL) {
  require("lmtest")
  require("sandwich")
  temp <- coeftest(x = lm_object, vcov. = vcovHC(lm_object, type = v_c_matrix))
  if (!is.null(var) & !is.null(var) & !is.null(group)) {
    return(make.df(var = var, comp = comp , n = length(lm_object$fitted.values),
                   model = temp, group = group))
  } else {
    return(temp)
  }
}

# plot function
make_plots <- function(res = res, group = "main", 
                       color_val = rev(c("red", "darkgreen", "blue")),
                       ncol_num = 2, make_facets = TRUE, xlab, ylab, 
                       comp_type = unique(res$comp), ggtitle = NULL) {
  temp <- ggplot(res[res$group %in% group & res$comp %in% comp_type,], 
                          aes(x=comp, y=coeff,color=comp,shape=comp)) +    
    geom_hline(yintercept = 0, colour = gray(0.5), lty = 2) +
    geom_linerange(aes(x = comp, ymin = coeff + qnorm(0.005)*se, 
                       ymax = coeff + qnorm(1 - 0.005)*se),
                   lwd = 1,position = position_dodge(0.5)) +
    geom_linerange(aes(x = comp, ymin = coeff + qnorm(0.025)*se, 
                       ymax = coeff + qnorm(1-0.025)*se),
                   lwd = 1.5,position = position_dodge(0.5)) +
    geom_point(fill = "white", size = 5) + 
    scale_colour_manual(values = color_val) + 
    scale_shape_manual(values = 21:23) +
    xlab(xlab) + ylab(ylab) + 
    coord_flip() + 
    theme_bw() + ggtitle(ggtitle) + 
    theme(legend.position="none") +
    theme(panel.grid.minor.x= element_line(colour=gray(0.9)),
          panel.grid.major.x= element_line(colour=gray(0.9))) + 
    theme(panel.grid.major.y = element_blank()) + 
    theme(axis.ticks.y = element_blank()) +
    theme(axis.text.y = element_text(face="bold")) +
    theme(plot.title=element_text(face="bold")) +
    theme(strip.text.x = element_text(face = "bold"))
  if (make_facets) {
    temp + facet_wrap(~var, ncol = ncol_num)  
  } else {
    temp
  }
}

# fix some variables
# mean impute CP$MagStake values where RepStake is missing
CP$MagStake[is.na(CP$MagStake) & CP$RepStake == 1] <- 
  round(mean(CP$MagStake, na.rm = TRUE))
# create new variable for reputation stakes that combines RepStake and MagStake
CP$RepStake_num <- NA
CP$RepStake_num[is.na(CP$MagStake)] <- 0
CP$RepStake_num[CP$MagStake == 1] <- 1
CP$RepStake_num[CP$MagStake == 2] <- 2
CP$RepStake_num[CP$MagStake == 3] <- 3

# standardize military intention variable
CP$MilInt_std <- CP$MilInt
CP$MilInt_std[CP$Incident.f %in% c("Attack", "Control")] <-
  scale(CP$MilInt[CP$Incident.f %in% c("Attack", "Control")])
CP$MilInt_std[CP$Incident.f %in% c("Collision", "Non-Collision")] <-
  scale(CP$MilInt[CP$Incident.f %in% c("Collision", "Non-Collision")])

# generate density plot 
CP$type <- NA
CP$type[CP$Incident.f %in% c("Attack", "Control")] <- "Basic"
CP$type[CP$Incident.f %in% c("Collision", "Non-Collision")] <- "ENE"
CP$treatment <- NA
CP$treatment[CP$Incident.f %in% c("Attack", "Collision")] <- "U.S. Pilot Fatality"
CP$treatment[CP$Incident.f %in% c("Control", "Non-Collision")] <- "No U.S. Pilot Fatality"
# get the group means
group_mean <- CP %>% group_by(type, treatment) %>% dplyr::summarise(
  gmean_MilInt = mean(MilInt)
)
# merge in with the data
CP <- merge(x = CP, y = group_mean, all.x = TRUE)
CP$treatment <- factor(CP$treatment, levels = c("U.S. Pilot Fatality", "No U.S. Pilot Fatality"))
ggplot(data = CP[!is.na(CP$type),], aes(x = MilInt, fill = treatment)) +
  geom_density(alpha = 0.5, adjust = 1.75) +
  geom_vline(mapping = aes(xintercept = gmean_MilInt, color = treatment),
             linetype = "longdash")  +
  facet_wrap(~type, ncol = 1) + 
  scale_x_continuous(labels = c("Very Unlikely", "Unlikely", "As Likely As Not",
                                "Likely", "Very Likely")) +
  theme_bb() + xlab("Likelihood of China Planned to Expand Presence and Capabilities") +
  scale_color_manual(values = c("orange3", "steelblue4"), name = "Treatment Assignment") +
  scale_fill_manual(values = c("orange1", "steelblue1"), name = "Treatment Assignment") +
  ylab("Density")
ggsave(paste0(images_directory, "Figure_33.pdf"), height=4.5, width=7.25,
       dpi = 300, colormodel = print_colormodel)
print("Complete: density plot")
Sys.time()

# Make regressions
res <- rbind(coef_r(lm_object = lm(MilInt ~ treatment == "U.S. Pilot Fatality", data=CP[CP$type == "Basic",]),
              var = "Hostile Military Intent", comp = "Basic", group = "IE"),
      coef_r(lm_object = lm(MilInt ~ treatment == "U.S. Pilot Fatality", data=CP[CP$type == "ENE",]),
             var = "Hostile Military Intent", comp = "ENE", group = "IE"),
      coef_r(lm_object = lm(MilInt_std ~ treatment == "U.S. Pilot Fatality", data=CP[CP$type == "Basic",]),
             var = "Hostile Military Intent", comp = "Basic", group = "IE_std"),
      coef_r(lm_object = lm(MilInt_std ~ treatment == "U.S. Pilot Fatality", data=CP[CP$type == "ENE",]),
             var = "Hostile Military Intent", comp = "ENE", group = "IE_std"))
      

res$comp <- factor(res$comp, levels = rev(levels(res$comp)))
make_plots(res = res, group = c("IE"), ncol_num = 1, color_val = c("red", "darkgreen"),
           ylab = "Difference-in-Means on a 5-Point Scale", xlab = NULL) 
ggsave(paste0(images_directory, "Figure_34.pdf"), height=2.5, width=5,
       dpi = 300, colormodel = print_colormodel)
make_plots(res = res, group = c("IE_std"), ncol_num = 1, color_val = c("red", "darkgreen"),
           xlab = NULL, ylab = "Standardized Difference-in-Means") 
ggsave(paste0(images_directory, "Figure_04.pdf"), height=2.5, width=5,
       dpi = 300, colormodel = print_colormodel)
print("Complete: standardized outcome coefficient plot")
Sys.time()

                       
